package cn.car91.utils;


public class EarthUtil {
    private final static double KmPerDegree=111.12000071117;
    private final static double DegreesPerKm=1.0/KmPerDegree;
    private final static double RadiansPerDegree=Math.PI/180.0;
    private final static double HalfPI =Math.PI/2;
    private final static double TwoPI =Math.PI*2;
    private final static double DegreesPerRadian=180.0/Math.PI;

    public static class Location {
        public double lat;
        public double lon;
    }
    
    /**
     * 
     * @param lat
     * @param lon
     * @param distance 单位公里
     * @return
     */
    public static Location getLocationRadius(double lat, double lon, double distance/*,double angle*/) {
        Location location = new Location();
        
        double  dist = distance;
        /*double  course = angle;*/           /* -> initial great-circle course (degrees) */
        double  course = 45;           /* -> initial great-circle course (degrees) */
        double  slt = lat;    /* -> starting decimal latitude (-S) */
        double  slg = lon;   /* -> starting decimal longitude(-W) */
        double  xlt = 0;    /* <- ending decimal latitude (-S) */
        double  xlg = 0;    /* <- ending decimal longitude(-W) */
        
        double  c, d, dLo, L1, L2, coL1, coL2, l;
        
        if (dist > KmPerDegree*180.0) {
            course -= 180.0;
            if (course < 0.0) course += 360.0;
            dist = KmPerDegree*360.0-dist;
        }
        
        if (course > 180.0) course -= 360.0;
        c = course*RadiansPerDegree;
        d = dist*DegreesPerKm*RadiansPerDegree;
        L1 = slt*RadiansPerDegree;
        slg *= RadiansPerDegree;
        coL1 = (90.0-slt)*RadiansPerDegree;
        
        coL2 = ahav(hav(c)/(1/(Math.cos(L1)*Math.sin(d)))+hav(d-coL1));
        L2 = HalfPI-coL2;
        l = L2-L1;
        if ((dLo=(Math.cos(L1)*Math.cos(L2))) != 0.0)
            dLo  = ahav((hav(d)-hav(l))/dLo);
        if (c < 0.0) dLo = -dLo;
        slg += dLo;
        if (slg < -Math.PI)
            slg += TwoPI;
        else if (slg > Math.PI)
            slg -= TwoPI;
        
        xlt = L2*DegreesPerRadian;
        xlg = slg*DegreesPerRadian;
        
        location.lat = Math.abs(lat - xlt);
        location.lon = Math.abs(lon - xlg);
        
        return location;
    }
    
    private static double hav(double x) {
        return ((1.0-Math.cos(x))*0.5);
    }
    
    private static double ahav(double x) {
        double t1,t2,t3;
        t1=1.0-(x*2.0);
        t2=(t1<0.0)?-1:1;
        t3=(Math.abs(t1)>1.0)?t2:t1;
        return Math.acos(t3);
    }
}
